<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Chebychev design of an FIR filter given a desired H(w)</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/filter_design/html/fir_chebychev_design.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Chebychev design of an FIR filter given a desired H(w)</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs an FIR filter given a desired frequency response H_des(w).</span>
<span class="comment">% The design is judged by the maximum absolute error (Chebychev norm).</span>
<span class="comment">% This is a convex problem (after sampling it can be formulated as an SOCP).</span>
<span class="comment">%</span>
<span class="comment">%   minimize   max |H(w) - H_des(w)|     for w in [0,pi]</span>
<span class="comment">%</span>
<span class="comment">% where H is the frequency response function and variable is h</span>
<span class="comment">% (the filter impulse response).</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% problem specs</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% number of FIR coefficients (including the zeroth one)</span>
n = 20;

<span class="comment">% rule-of-thumb frequency discretization (Cheney's Approx. Theory book)</span>
m = 15*n;
w = linspace(0,pi,m)'; <span class="comment">% omega</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% construct the desired filter</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% fractional delay</span>
D = 8.25;            <span class="comment">% delay value</span>
Hdes = exp(-j*D*w);  <span class="comment">% desired frequency response</span>

<span class="comment">% Gaussian filter with linear phase (uncomment lines below for this design)</span>
<span class="comment">% var = 0.05;</span>
<span class="comment">% Hdes = 1/(sqrt(2*pi*var))*exp(-(w-pi/2).^2/(2*var));</span>
<span class="comment">% Hdes = Hdes.*exp(-j*n/2*w);</span>

<span class="comment">%*********************************************************************</span>
<span class="comment">% solve the minimax (Chebychev) design problem</span>
<span class="comment">%*********************************************************************</span>
<span class="comment">% A is the matrix used to compute the frequency response</span>
<span class="comment">% A(w,:) = [1 exp(-j*w) exp(-j*2*w) ... exp(-j*n*w)]</span>
A = exp( -j*kron(w,[0:n-1]) );

<span class="comment">% optimal Chebyshev filter formulation</span>
cvx_begin
  variable <span class="string">h(n,1)</span>
  minimize( max( abs( A*h - Hdes ) ) )
cvx_end

<span class="comment">% check if problem was successfully solved</span>
disp([<span class="string">'Problem is '</span> cvx_status])
<span class="keyword">if</span> ~strfind(cvx_status,<span class="string">'Solved'</span>)
  h = [];
<span class="keyword">end</span>

<span class="comment">%*********************************************************************</span>
<span class="comment">% plotting routines</span>
<span class="comment">%*********************************************************************</span>
<span class="comment">% plot the FIR impulse reponse</span>
figure(1)
stem([0:n-1],h)
xlabel(<span class="string">'n'</span>)
ylabel(<span class="string">'h(n)'</span>)

<span class="comment">% plot the frequency response</span>
H = [exp(-j*kron(w,[0:n-1]))]*h;
figure(2)
<span class="comment">% magnitude</span>
subplot(2,1,1);
plot(w,20*log10(abs(H)),w,20*log10(abs(Hdes)),<span class="string">'--'</span>)
xlabel(<span class="string">'w'</span>)
ylabel(<span class="string">'mag H in dB'</span>)
axis([0 pi -30 10])
legend(<span class="string">'optimized'</span>,<span class="string">'desired'</span>,<span class="string">'Location'</span>,<span class="string">'SouthEast'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling sedumi: 1199 variables, 321 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 321, order n = 901, dim = 1200, blocks = 300
nnz(A) = 12580 + 0, nnz(ADA) = 13341, nnz(L) = 6831
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.03E+02 0.000
  1 :  -1.78E+00 9.13E+01 0.000 0.4496 0.9000 0.9000   1.40  1  1  2.6E+02
  2 :  -8.52E-01 3.66E+01 0.000 0.4009 0.9000 0.9000   4.03  1  1  3.4E+01
  3 :  -6.50E-01 2.07E+01 0.000 0.5648 0.9000 0.9000   2.79  1  1  1.3E+01
  4 :  -7.18E-01 6.34E+00 0.000 0.3065 0.9000 0.9000   1.19  1  1  4.2E+00
  5 :  -7.07E-01 3.37E-01 0.000 0.0532 0.9900 0.9900   1.11  1  1  2.1E-01
  6 :  -7.07E-01 1.72E-02 0.135 0.0510 0.9900 0.9900   1.00  1  1  1.1E-02
  7 :  -7.07E-01 2.56E-03 0.000 0.1491 0.9042 0.9000   1.00  1  1  1.6E-03
  8 :  -7.07E-01 6.97E-05 0.000 0.0272 0.9900 0.0000   1.00  1  1  5.4E-05
  9 :  -7.07E-01 6.98E-08 0.000 0.0010 0.9990 0.9990   1.00  1  1  6.5E-08
 10 :  -7.07E-01 6.96E-10 0.000 0.0100 0.9990 0.9939   1.00  1  2  6.6E-10

iter seconds digits       c*x               b*y
 10      0.1   Inf -7.0710678098e-01 -7.0710678095e-01
|Ax-b| =   4.8e-10, [Ay-c]_+ =   1.2E-10, |x|=  1.7e+00, |y|=  5.1e+00

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    8.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=2, |skip| = 0, ||L.L|| = 2591.85.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.707107
Problem is Solved
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_chebychev_design__01.png" alt=""> <img src="fir_chebychev_design__02.png" alt=""> 
</div>
</div>
</body>
</html>